À la fin de ce chapitre, vous
Que voulez-vous exprimer?
Que voulez-vous exprimer?
Disposition des trous de balle sur les avions britaniques de retour de leur sortie durant la deuxième guerre mondiale.
Source: Whitlock et Schuluter (2015), The Analysis of Biological Data.
La statistique est l’étude des méthodes pour mesurer des aspects de populations à partir d’échantillons et pour quantifier l’incertitude des mesures. (ma traduction de Whitlock et Schuluter (2015)).
Ensemble circonscrit du sujet d’étude.
Sous-ensemble de la population.
Si l’échantillon est représentatif de la population (non biaisé), on peut extrapoler ses propriétés à la population (inférence).
Dimension d’un objet catractérisé par une distribution de probabilité.
Une probabilité est la vraisemblance qu’un évènements se réalise chez un échantillon.
Une distribution décrit la probabilité d’obtenir une valeur (distribution discrète) ou une plage de valeurs (distribution continue) dans une échantillon pris au hasard d’une population.
Une distribution des probabilité est décrite par un ou plusieurs paramètres, par exemple une moyenne et un écart-type dans le cas d’une distribution normale.
Une statistique est une estimation d’un paramètre calculée à partir des données, par exemple une moyenne et un écart-type échantillonnaux.
Fréquentielle. Les données sont générées par des mécanismes stochastiques décrits par des distributions de probabilités, et nous cherchons à déterminer la probabilité que les données soient générées par ces mécanismes (p-valule). <– approche la plus commune (et couverte dan ce chapitre)
Bayésienne. Les paramètres et leur incertitude sont déterminés par les données et les connaissances préalables. <– approche de plus en plus utilisée (effleurée dans le chapitre 6, en extra)
Tout dépend de la question posée.
exemple tirée du blogue Dynamic Ecology
Un test d’hypothèse permet de décider si une hypothèse est confirmée ou rejetée à un seuil de probabilité prédéterminé.
H0. L’hypothèse nulle propose l’absence d’effet statistique (c’est l’hypothèse de l’avocat du diable 😈) .
H1. L’hypothèse alternative propose l’inverse (😇).
Est-ce que les données sont conforment avec avec l’hypothèse de la vie sur Mars.
H0. Il n’y a pas de vie sur Mars.
H1. Il y a de la vie sur Mars.
Dans une expérience, la confirmation de l’hypothèse nulle n’est pas un échec.
Comparaison des moyennes deux échantillons dont la distribution est normale et dont la variance est la même.
## # A tibble: 5 x 2
## echantillon value
## <chr> <dbl>
## 1 a 10.1
## 2 b 7.98
## 3 a 5.82
## 4 b 11.8
## 5 a 9.55
tt_exemple <- t.test(formula = value ~ echantillon,
data = exemple_ech, var.equal = TRUE)
tt_exemple##
## Two Sample t-test
##
## data: value by echantillon
## t = 1.5581, df = 58, p-value = 0.1247
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2990577 2.3995879
## sample estimates:
## mean in group a mean in group b
## 9.296136 8.245871
## List of 9
## $ statistic : Named num 1.56
## ..- attr(*, "names")= chr "t"
## $ parameter : Named num 58
## ..- attr(*, "names")= chr "df"
## $ p.value : num 0.125
## $ conf.int : num [1:2] -0.299 2.4
## ..- attr(*, "conf.level")= num 0.95
## $ estimate : Named num [1:2] 9.3 8.25
## ..- attr(*, "names")= chr [1:2] "mean in group a" "mean in group b"
## $ null.value : Named num 0
## ..- attr(*, "names")= chr "difference in means"
## $ alternative: chr "two.sided"
## $ method : chr " Two Sample t-test"
## $ data.name : chr "value by echantillon"
## - attr(*, "class")= chr "htest"
## [1] 0.1246579
“La p-value n’a jamais été conçue comme substitut au raisonnement scientifique” Ron Wasserstein, directeur de l’American Statistical Association [ma traduction].
Un résultat montrant une p-value plus élevée que 0.05 est-il pertinent?
Lors d’une conférence, Dr Evil ne présentent que les résultats significatifs de ses essais au seuil de 0.05. Certains essais ne sont pas significatifs, mais bon, ceux-ci ne sont pas importants…
Dans une expérience, la confirmation de l’hypothèse nulle n’est pas un échec.
Modifier les données ou la question posée dans le but d’obtenir une p-value < 0.05, c’est tricher.
Comparaison des moyennes de plus de 2 groupes (à variance égale)
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 437.1 218.55 1180 <2e-16 ***
## Residuals 147 27.2 0.19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\[Y = f\left( X \right)\] \(X\): variables prédictives / indépendantes / covariables / explicatives \(Y\): variables de sortie / réponse / dépendante
| ⚠️ Ne pas confondre. ⚠ |
Modèles prédictifs. Le modèle doit être testé sur des données qui n’ont pas servi à l’entraîner. (chapitre 12)
Modèle explicatif. Le modèle sert à tester des hypothèses. (chapitres 5)
Variables fixes: testées lors de l’expérience: dose du traitement, espèce/cultivar, météo, etc.
Variables aléatoires: sources de variation qui génèrent du bruit dans le modèle: les unités expérimentales ou le temps lors de mesures répétées.
\[y = \beta_0 + \beta_1 x + \epsilon\]
Question. Le rendement varie-t-il selon la dose d’azote?
H0. la dose d’azote n’affecte pas le rendement (\(\beta_1 = 0\)).
lm##
## Call:
## lm(formula = yield ~ nitro, data = lasrosas.corn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -53.183 -15.341 -3.079 13.725 45.897
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 65.843213 0.608573 108.193 < 2e-16 ***
## nitro 0.061717 0.007868 7.845 5.75e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.66 on 3441 degrees of freedom
## Multiple R-squared: 0.01757, Adjusted R-squared: 0.01728
## F-statistic: 61.54 on 1 and 3441 DF, p-value: 5.754e-15
predictggplot(data = lasrosas.corn, mapping = aes(x = nitro, y = yield)) +
geom_point() +
geom_line(data = tibble(nitro = lasrosas.corn$nitro, yield = predict(modlin_1)),
size = 2, colour = "grey35")Résidus: erreurs du modèle, \(\epsilon = y - \hat{y}\), fonction residuals().
res_df <- data.frame(nitro = lasrosas.corn$nitro,
residus_lm = residuals(modlin_1))
res_df %>% sample_n(4)## nitro residus_lm
## 3086 99.8 -14.27259
## 231 106.0 -16.50523
## 2911 124.6 -35.33317
## 800 106.0 -18.80523
\[ y = X \beta + \epsilon \]
topo est la classe topographique
## year lat long yield nitro topo bv rep nf
## 1 1999 -33.05113 -63.84886 72.14 131.5 W 162.60 R1 N5
## 35 1999 -33.05155 -63.84645 56.70 131.5 HT 169.49 R1 N5
## 51 1999 -33.05174 -63.84532 59.94 131.5 E 182.12 R1 N5
## 67 1999 -33.05195 -63.84412 70.35 131.5 LO 168.00 R1 N5
## (Intercept) topoHT topoLO topoW
## 1 1 0 0 1
## 35 1 1 0 0
## 51 1 0 0 0
## 67 1 0 1 0
## topoE topoHT topoLO topoW
## 1 0 0 0 1
## 35 0 1 0 0
## 51 1 0 0 0
## 67 0 0 1 0
##
## Call:
## lm(formula = yield ~ topo, data = lasrosas.corn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.371 -11.933 -1.593 11.080 44.119
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 78.6653 0.5399 145.707 <2e-16 ***
## topoHT -30.0526 0.7500 -40.069 <2e-16 ***
## topoLO 6.2832 0.7293 8.615 <2e-16 ***
## topoW -11.8841 0.7039 -16.883 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.59 on 3439 degrees of freedom
## Multiple R-squared: 0.4596, Adjusted R-squared: 0.4591
## F-statistic: 975 on 3 and 3439 DF, p-value: < 2.2e-16
## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## topo 3 622351 207450 974.96 < 2.2e-16 ***
## Residuals 3439 731746 213
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Intercept) statut_o.L statut_o.Q statut_o.C statut_o^4
## 1 1 -0.6324555 0.5345225 -3.162278e-01 0.1195229
## 2 1 -0.3162278 -0.2672612 6.324555e-01 -0.4780914
## 3 1 0.0000000 -0.5345225 -4.095972e-16 0.7171372
## 4 1 0.3162278 -0.2672612 -6.324555e-01 -0.4780914
## 5 1 0.6324555 0.5345225 3.162278e-01 0.1195229
## attr(,"assign")
## [1] 0 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$statut_o
## [1] "contr.poly"
modlin_mult <- lm(yield ~ lat + long + nitro + topo + bv,
data = lasrosas.corn)
summary(modlin_mult)##
## Call:
## lm(formula = yield ~ lat + long + nitro + topo + bv, data = lasrosas.corn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.405 -11.071 -1.251 10.592 40.078
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.946e+05 3.309e+04 5.882 4.45e-09 ***
## lat 5.541e+03 4.555e+02 12.163 < 2e-16 ***
## long 1.776e+02 4.491e+02 0.395 0.693
## nitro 6.867e-02 5.451e-03 12.597 < 2e-16 ***
## topoHT -2.665e+01 1.087e+00 -24.520 < 2e-16 ***
## topoLO 5.565e+00 1.035e+00 5.378 8.03e-08 ***
## topoW -1.465e+01 1.655e+00 -8.849 < 2e-16 ***
## bv -5.089e-01 3.069e-02 -16.578 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.47 on 3435 degrees of freedom
## Multiple R-squared: 0.5397, Adjusted R-squared: 0.5387
## F-statistic: 575.3 on 7 and 3435 DF, p-value: < 2.2e-16
modlin_mult_sc <- lm(yield ~ lat + long + nitro + topo + bv,
data = lasrosas.corn_sc)
summary(modlin_mult_sc)##
## Call:
## lm(formula = yield ~ lat + long + nitro + topo + bv, data = lasrosas.corn_sc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.405 -11.071 -1.251 10.592 40.078
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 78.9114 0.6666 118.376 < 2e-16 ***
## lat 3.9201 0.3223 12.163 < 2e-16 ***
## long 0.3479 0.8796 0.395 0.693
## nitro 2.9252 0.2322 12.597 < 2e-16 ***
## topoHT -26.6487 1.0868 -24.520 < 2e-16 ***
## topoLO 5.5647 1.0347 5.378 8.03e-08 ***
## topoW -14.6487 1.6555 -8.849 < 2e-16 ***
## bv -4.9253 0.2971 -16.578 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.47 on 3435 degrees of freedom
## Multiple R-squared: 0.5397, Adjusted R-squared: 0.5387
## F-statistic: 575.3 on 7 and 3435 DF, p-value: < 2.2e-16
Une intéraction est une pente additionnelle informant sur l’effet statistique combiné de deux variables.
Interface-formule: le symbole : appelle l’intéraction et le synbole * appelle les effets simples et les intéractions.
##
## Call:
## lm(formula = yield ~ nitro * topo, data = lasrosas.corn_sc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.984 -11.985 -1.388 10.339 40.636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 78.6999 0.5322 147.870 < 2e-16 ***
## nitro 1.8131 0.5351 3.388 0.000711 ***
## topoHT -30.0052 0.7394 -40.578 < 2e-16 ***
## topoLO 6.2026 0.7190 8.627 < 2e-16 ***
## topoW -11.9628 0.6939 -17.240 < 2e-16 ***
## nitro:topoHT 1.2553 0.7461 1.682 0.092565 .
## nitro:topoLO 0.5695 0.7186 0.792 0.428141
## nitro:topoW 0.7702 0.6944 1.109 0.267460
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.38 on 3435 degrees of freedom
## Multiple R-squared: 0.4756, Adjusted R-squared: 0.4746
## F-statistic: 445.1 on 7 and 3435 DF, p-value: < 2.2e-16
Évaluer l’effet sur une variable entière, booléenne, factorielle, etc.
glm##
## Call:
## glm(formula = worms ~ trt, family = "poisson", data = cochran.wireworms)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8279 -0.9455 -0.2862 0.6916 3.1888
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1823 0.4082 0.447 0.655160
## trtM 1.6422 0.4460 3.682 0.000231 ***
## trtN 1.7636 0.4418 3.991 6.57e-05 ***
## trtO 1.5755 0.4485 3.513 0.000443 ***
## trtP 1.3437 0.4584 2.931 0.003375 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 64.555 on 24 degrees of freedom
## Residual deviance: 38.026 on 20 degrees of freedom
## AIC: 125.64
##
## Number of Fisher Scoring iterations: 5
\[ y = A \left( 1 - e^{-R \left( E + x \right)} \right) \]
nlsmodnl_1 <- nls(yield ~ A * (1 - exp(-R*(E + nitro))),
data = engelstad.nitro,
start = list(A = 75, E = 30, R = 0.02))
summary(modnl_1)##
## Formula: yield ~ A * (1 - exp(-R * (E + nitro)))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A 75.023427 3.331860 22.517 <2e-16 ***
## E 66.164111 27.251592 2.428 0.0184 *
## R 0.012565 0.004881 2.574 0.0127 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.34 on 57 degrees of freedom
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 8.063e-06
À la différence d’un effet fixe, un effet aléatoire (\(b\)) sur des variables aléatoires (\(Z\)) sera toujours distribué normalement avec une moyenne de 0 et une certaine variance. Pour un modèle à intercept aléatoire, nous aurons
\[ y = X \beta + Z b + \epsilon \]
Pour un modèle à pente aléatoire, nous aurons (à vérifier)
\[ y = X \left( \beta + Zb \right) + \epsilon \]
lme du module nlmerandom = ~ slope|intercept
mmodlin_1 <- lme(fixed = yield ~ lat + long + nitro + topo + bv,
random = ~ 1|year/rep,
data = lasrosas.corn)
summary(mmodlin_1)## Linear mixed-effects model fit by REML
## Data: lasrosas.corn
## AIC BIC logLik
## 26535.37 26602.93 -13256.69
##
## Random effects:
## Formula: ~1 | year
## (Intercept)
## StdDev: 20.35425
##
## Formula: ~1 | rep %in% year
## (Intercept) Residual
## StdDev: 11.17447 11.35617
##
## Fixed effects: yield ~ lat + long + nitro + topo + bv
## Value Std.Error DF t-value p-value
## (Intercept) -1379436.9 55894.55 3430 -24.679273 0.000
## lat -25453.0 1016.53 3430 -25.039084 0.000
## long -8432.3 466.05 3430 -18.092988 0.000
## nitro 0.0 0.00 3430 1.739757 0.082
## topoHT -27.7 0.92 3430 -30.122438 0.000
## topoLO 6.8 0.88 3430 7.804733 0.000
## topoW -16.7 1.40 3430 -11.944793 0.000
## bv -0.5 0.03 3430 -19.242424 0.000
## Correlation:
## (Intr) lat long nitro topoHT topoLO topoW
## lat 0.897
## long 0.866 0.555
## nitro 0.366 0.391 0.247
## topoHT 0.300 -0.017 0.582 0.024
## topoLO -0.334 -0.006 -0.621 -0.038 -0.358
## topoW 0.403 -0.004 0.762 0.027 0.802 -0.545
## bv -0.121 -0.012 -0.214 -0.023 -0.467 0.346 -0.266
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.32360269 -0.66781575 -0.07450856 0.61587533 3.96434001
##
## Number of Observations: 3443
## Number of Groups:
## year rep %in% year
## 2 6
## Level: year
## (Intercept)
## 1999 -13.71488
## 2001 13.71488
##
## Level: rep %in% year
## (Intercept)
## 1999/R1 -12.3828119
## 1999/R2 -2.0608423
## 1999/R3 10.3099880
## 2001/R1 -9.2956080
## 2001/R2 0.8470003
## 2001/R3 12.5822738
## mean(`(Intercept)`)
## 1 -2.279475e-11
mm <- nlme(yield ~ A * (1 - exp(-R*(E + nitro))),
data = engelstad.nitro,
start = c(A = 75, E = 30, R = 0.02),
fixed = list(A ~ 1, E ~ 1, R ~ 1),
random = A ~ year | loc)
summary(mm)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: yield ~ A * (1 - exp(-R * (E + nitro)))
## Data: engelstad.nitro
## AIC BIC logLik
## 477.2285 491.8889 -231.6143
##
## Random effects:
## Formula: A ~ year | loc
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## A.(Intercept) 2.602881636 A.(In)
## A.year 0.003066058 -0.556
## Residual 11.152760033
##
## Fixed effects: list(A ~ 1, E ~ 1, R ~ 1)
## Value Std.Error DF t-value p-value
## A.(Intercept) 74.58277 4.722806 56 15.792046 0.0000
## E 65.57006 25.534747 56 2.567876 0.0129
## R 0.01308 0.004807 56 2.720245 0.0087
## Correlation:
## A.(In) E
## E 0.379
## R -0.483 -0.934
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.83376801 -0.89290179 0.07419027 0.68354461 1.82431474
##
## Number of Observations: 60
## Number of Groups: 2
À la fin de ce chapitre, vous